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Abstract 

Incomplete lineage sorting (ILS) is a common source of gene tree in- 
congruence in multilocus analyses. A large number of methods have been 
developed to infer species trees in the presence of ILS. Here we provide a 
mathematical analysis of several coalescent-based methods. Our analysis is 
performed on a three-taxon species tree and assumes that the gene trees are 
correctly reconstructed along with their branch lengths. 

1 Introduction 

Incomplete lineage sorting (ILS) is an important confounding factor in phyloge- 
netic analyses based on multiple genes or loci [ Mad97. DR09] . ILS is a population- 
level phenomenon that is caused by the failure of two lineages to coalesce in a 
population, leading to the possibility that one of the lineages first coalesces with 
a lineage from a less closely related population. As a result, it can produce exten- 
sive gene tree incongruence that must be accounted for appropriately in multilocus 
analyses IIDR06L 

A large number of methods have been developed to address this source of in- 
congruence llLYK+0911 . Several such methods rely on a statistical model of ILS 
known as the multispecies coalescent. In this model, populations are connected by 
a phylogeny. Independent coalescent processes are performed in each population 
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and assembled to produce gene trees. Several methods have been shown to be sta- 
tistically consistent under the multispecies coalescent, that is, they are guaranteed 
to return the correct species tree given enough loci. 

The performance and accuracy of coalescent-based multilocus methods have 
been the subject of numerous simulation studies ULRlli |LYPE09[ lYNlOL In 
this paper, we complement such studies with a detailed analytical comparison 
in a tractable test case, a three-taxon species tree. We analyze 7 methods: max- 
imum likelihood (ML), GLASS/Maximum Tree (MT), R*, STAR, minimizing 
deep coalescences (MDC), STEAC, and shallowest coalescences (SC). Under the 
assumption that gene trees are reconstructed without estimation error, we derive 
the exponential decay rate of the failure probability as the internal branch length 
of the species tree varies. The analysis, which relies on large-deviations theory, 
reveals that ML and GLASS/MT are more accurate in this setting than the other 
methods — especially in the regime where ILS is more common. 

2 Materials and Methods 

2.1 Multispecies coalescent: Three-taxon case 

We first describe the statistical model under which our analysis is performed, the 
multispecies coalescent. We only discuss the three-taxon case. For more details, 
see [IDR09II and references therein. 

A weighted rooted tree is called ultrametric if each leaf is exactly at the same 
distance from the root. For a three-leaf ultrametric tree G with leaves a, b, and c, 
we denote by ab\c the topology where a and b are closer to each other than to c, 
and similarly for ac\b, bc\a. The topology of G is denoted by T[G]. 

Let S be an ultrametric species phylogeny with three taxa. We assume that 
all haploid populations in S have population size N. We denote the current pop- 
ulations by A, B and C (which we identify with the leaves of S) and we assume 
that S has topology AB|C. The ancestral populations are AB (corresponding 
to the immediate ancestor to populations A and B) and ABC (corresponding to 
the ancestor of populations A, B and C). The corresponding divergence times 
(backwards in time from the present) are denoted by tab and tabc with the as- 
sumption Tab < Tabc- All times are given in units of generations. For a 
population X, we let be the divergence time of the parent population of X. Let 
X = {A, B, C, AB, ABC} be the set of all populations in S. 

We consider L loci i = 1, . . . , L and, for each locus, we sample one lineage 
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from each population at time 0. For locus we denote by the number of 
lineages entering population X and by O^^ the number of lineages exiting pop- 
ulation X (backwards in time), where necessarily I^^ > Similarly, for 
k = O^^ + 1, . . . , the time of the coalescent event bringing the number of 
lineages from to A; — 1 in population X and locus i is T^'''\ We denote by 
Gi, . . . , Gl the corresponding ultrametric gene trees (including both topology and 
branch lengths). 

Then, under the multispecies coalescent, assuming the loci are unlinked, the 
likelihood of the gene trees is given by 

e=i ^ xex V / \ / 

(il^^^'+l) H h 

where we let T-^' ^ = tx for convenience nRY03ll . 

The parameter governing the extent of incomplete lineage sorting is the length 
of the internal branch of S 

t = tabc — Tab- 

The probability that the lineages from A and B fail to coalesce in branch AB, an 
event we denote by Fail^ for locus i (and its complement by SUCCESS^), is 

1 - p = e"*. 

Note that, in that case, all three gene-tree topologies are equally likely. Of course, 
1-p^-last-^O. 

2.2 Multilocus methods 

A basic goal of multilocus analyses is to reconstruct a species phylogeny (includ- 
ing possibly estimates of the divergence times) from a collection of gene trees. 
Here we assume that the data consists of L gene trees Gi, . . . ,Gl corresponding 
to L unlinked loci generated under the multispecies coalescent. We assume fur- 
ther that the gene trees are ultrametric and that their topology and branch lengths 
are estimated without error. 
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We consider several common multilocus methods. In our setting, several of 
these methods are in fact equivalent and we therefore group them below. Note 
further that we only consider statistically consistent methods, that is, methods that 
are guaranteed to converge on the right species phylogeny as the number of loci L 
increases to +00 (at least, in the test case we described above). We briefly describe 
these methods. For more details, see e.g. [lLYK+091 and references therein. 



ML/GLASS/MT Under the multispecies coalescent, maximum likelihood (ML) 
selects the topology and divergence times that maximizes the likelihood ([T]). 

In the GLASS method P MRIOI . the species phylogeny is reconstructed from 
a distance matrix in which the entries are the minimum gene coalescence times 
across loci. The equivalent Maximum Tree (MT) method was introduced and 
studied in [iLPOTllErPOTllLYPlOII . 

A key result in ULYPIOH is that, in the constant-population case, the term inside 
the exponential in the likelihood ([T]) is monotonically decreasing in the divergence 
times. As a result, because GLASS and MT select the phylogeny with the largest 
possible divergence times, maximum likelihood is equivalent to GLASS and MT 
in this context. See ULYPIOH for details. 



i?*/STAR/MDC In the R* consensus method gBry03[lDDBR09ll . for each three- 
taxon set (here, we only have one such set), we include the topology that appears in 
highest frequency among the loci and we reconstruct the most resolved phylogeny 
that is compatible with these three-taxon topologies. 

In the STAR method [ILYPE09II , the species phylogeny is reconstructed from 
a distance matrix in which the entries are the average ranks of gene coalescence 
times across loci. Here the root has the highest rank and the rank decreases by one 
as one goes from the root to the leaves. 

The minimizing deep coalescences (MDC) method [|Mad97[|TN09l selects the 
species phylogeny that requires the smallest number of "extra lineages," that is, 
lineages that fail to coalesce in a branch of the species phylogeny. 

On a three-taxon phylogeny, there is only three distinct rooted topologies. In 
each case, the most recent divergence is assigned rank 1 in STAR and the other 
divergence is assigned rank 2. Hence selecting the topology corresponding to the 
lowest average rank is equivalent to selecting the most common topology among 
all loci — which is what R* does. A similar argument shows that MDC also selects 
the R* consensus tree in our test case. 
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STEAC/SC In the STEAC method [ILYPE09II . the species phylogeny is recon- 
structed from a distance matrix in which the entries are the average coalescence 
times across loci. The shallowest coalescences (SC) method is similar to STEAC 
in that it uses average coalescence times. The difference between the two methods 
is in how they deal with multiple alleles per population. Since we only consider 
the single-allele case, the two methods are equivalent here. 



2.3 Large-deviations approach 

As mentioned above, we consider estimation methods that are statistically con- 
sistent in the sense that they are guaranteed to converge on the correct species 
phylogeny as the number of loci L increases to +00. To compare different meth- 
ods, we derive the rate of exponential decay of the probability of failure. Let S be 
a species phylogeny with internal branch length t and assume that Gi, . . . , are 
unlinked gene trees generated under the multispecies coalescent. As L — > +00, 
large-deviations theory (see e.g. [|Dur96l ) gives a characterization of the (expo- 
nential) decay rate 

1 

amit) = - lim — In P [Method M fails given L loci from S] . 

That is, roughly 

P[Method M fails given L loci from S] ^ e"-^°"(*\ 

for large L. As the notation indicates, the key parameter that influences the decay 
rate is the length of the internal branch t of the species phylogeny. In particular, we 
expect that amit) is increasing in t as a larger t makes the reconstruction problem 
easier. 

To derive «m(^)5 we express the probability of failure as a large deviation event 
of the form 



P [Method M fails given L loci from S]=F 



where y is a constant and {Yi}f^^ are independent identically distributed random 
variables. The particular choice of random variables depends on the method, as 
we explain below. Let 



is) = E[e^^^" 
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be the moment-generating function of (which does not depend on i by assump- 
tion). Then the decay rate is given by 

amit) = ys^ - ln0(s*), (2) 

where s* > is the solution (if it exists) to 

provided there is an s > such that 0(s) < +oo, y > E[Y£] and Yi is not a point 
mass at E[Ye]. For more details on large-deviations theory, see e.g. [|Dur96L 

3 Results 

3.1 A domination result 

We first argue that, given perfectly reconstructed unhnked gene trees under the 
multispecies coalescent, ML/GLASS/MT always has a greater probability of suc- 
cess than i?*/STAR/MDC and STEAC/SC— or, in fact, any other method. Indeed 
note that the probability of success can be divided into two cases: 

1. The case where SUCCESS^ occurs for at least one locus i, an event of prob- 
ability (1 - (1 - p)^). In that case, ML/GLASS/MT necessarily succeeds 
whereas the other two methods succeed with probability < 1. 

2. The case where Fail^ occurs for all loci i, an event of probability (1 — p)^. 
In that case, all methods succeed with probability 1/3 by symmetry. For 
instance, for ML/GLASS/MT, any pair of populations is equally likely to 
lead to the smallest inter-species distance. A similar argument applies to 
the other two methods. 

Hence, overall ML/GLASS/MT succeeds with greater probabihty. 

3.2 Decay rates 

We derive the decay rates for the methods above. The results are plotted in Fig- 
ure [T] The asymptotic regimes are highlighted in Figures |2] and [3] All proofs can 
be found in the appendix. 
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ML/GLASS/MT In this case, the decay rate can be derived directly without 
using ([2]). Following the derivation in HMRIOI (see also ULYPIOII for a similar 
argument), ML/GLASS/MT succeeds with probability 

Then we get the following: 

Claim 1 (ML/GLASS/MT) The decay rate of MUGLASS/MT on S is 

i?*/STAR/MDC For a locus £, we let zfl be 1 if Fail^ occurs and r[Gi] = 
AB|C, and otherwise. We let 

L 

Similarly, we define zjj,, zl^^, Zac and Zbc- Then i?*/STAR/MDC fails if 

Zab + {L - Zac - ^bc - ^ab) < max{ZAc, ^bc}- 
It can be shown that 

aR*(t) = - lim ylnP[2ZAc + ^BC >^]- 

Then we get the following: 

Claim 2 (/?*/STAR/MDC) The decay rate of R*/STAR/MDC on S is 

Ast^Q, 

and, as t +00, 

/X t 1,4 
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STEAC/SC For a locus £, we let D^f^ be the time to the most recent common 
ancestor of A and B in Gc (in units of generations). We let 



Similarly, we define D^, D^^l, Vac and Vbc- Then STEAC/SC fails if 

Dab > min{DAc,^^Bc}- 

It can be shown that 

asTEAc(t) = lim lnP[pAB - ^^ac > 0]. 
Then we get the following: 

Claim 3 (STEAC/SC) The decay rate of STEAC/SC on S is 

"STEAcW = -In' 



3(1 -s2) 

where < s=k < 1 /5 the unique solution to the fixed-point equation 

= ^[6s.-3t(l-s^)]e(i--)*. 

Further, ast^O, 

asTEAcW = lt^ + 0{t^), 

o 

and, as t +00, 

asTEAc(^) ~ t — Int — 0.1656. 



4 Discussion 

As can be seen from Figures [T||2]and[3]as well as from the asymptotics, ML/GLA- 
SS/MT does indeed give a larger decay rate for all t. In fact, the decay rate of 
ML/GLASS/MT is significantly higher, especially as t — t- that is, under high 
levels of incomplete lineage sorting. For instance, to be concrete, if L = 500 loci 
and t = 0.1 (in units of N generations), the probability of failure is approximately: 
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1.9 X 10-22 for ML/GLASS/MT; 0.038 for i?*/STAR/MDC; 0.16 for STEAC/SC. 
Intuitively, this difference in behavior arises from the fact that ML/GLASS/MT re- 
quires only one successful locus, whereas i?*/STAR/MDC and STEAC/SC rely on 
an average over all loci. 

Comparing i?*/STAR/MDC and STEAC/SC, note that or. (t) is higher than 
asTEAc(^) for small t but that the situation is reversed for large t. In fact, in 
the limit t — )• +oo, asTEAc(^) grows at roughly the same rate as the optimal 
«ml(^)- At large t, STEAC/SC has somewhat of an advantage in that the expec- 
tation gap in the failure event increases linearly with t, whereas it saturates under 
i?*/STAR/MDC. 

The analysis described here ignores several features that influence the accu- 
racy of species tree reconstruction. Notably we have assumed that gene trees, 
including their branch lengths, are reconstructed without error. On real sequence 
datasets, the uncertainty arising from gene-tree estimation plays an important role. 
For instance, although GLASS/MT achieves the optimal decay rate in our setting, 
these methods are in fact sensitive to sequence noise because they rely on the 
computation of a minimum over loci — the very feature that leads to their superior 
performance here. Extending our analysis to incorporate gene tree estimation er- 
ror is an important open problem which should help in the design of multilocus 
methods. It is important to note that, under appropriate modeling of sequence data, 
ML is not in general equivalent to GLASS/MT and is likely to be more robust to 
estimation error. In particular our analysis suggest that ML may be significantly 
more accurate than other methods in multilocus studies. 

Other extensions deserve further study. Often many alleles are sampled from 
each population. Note that the benefit of multiple alleles is known to saturate 
as the number of alleles increases [|Ros02|| . This is because the probability of 
observing any number of alleles at the top of a branch is uniformly bounded in the 
number alleles existing at the bottom. 

Further, the molecular clock assumption, although it may be a reasonable first 
approximation in the context of recently diverged populations, should not be nec- 
essary for our analysis. One should also consider larger numbers of taxa, varying 
population sizes, etc. 

Simulation studies may provide further insight into these issues. However an 
analytical approach, such as the one we have used here, is valuable in that it allows 
the study of an entire class of models in one analysis. It can also provide useful, 
explicit predictions to guide the design of reconstruction procedures. 
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Figure 2: Decay rates as t — > 0. The dotted lines indicate the respective predicted 
asymptotics. 
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A Proofs 



A.l ML/GLASS/MT 

Decay rate Following the derivation in HMRIOII (see also ULYPIOII for a similar 
argument), ML/GLASS/MT succeeds with probability 



where the two terms correspond to the two cases described in Section 3.1 Hence 
the decay rate of the failure probability is 

aML(t) = lim -^ln(l-(l-(l-p)^) + J(l-p)^ 

1 f 2 r 

= lim - ln<^-(l-p)^ 

r 1 2 

= lim <— — In ln(l— p) 

L-s>+oo \^ L 3 

= — ln(l — p) 
= t. 



A.2 i?7STAR/MDC 

Definitions For a locus £, we let Z% be 1 if Fail^ occurs and T[G(\ = AB|C, 
and otherwise. We let 

L 

^AB = ^ ^AB- 

Similarly, we define zjj,, zjf^, Zac and Zbc- Then i?*/STAR/MDC fails if 

2ab + {L - Zac - ^bc - ^ab) < max{ZAc, ^bc}, 

an event we denote by £. The second term on the LHS comes form the fact that, 
given Success^, T[Gi\ = AB|C. To deal with the term on the RHS, we re-write 
£ as 

2 max{ZAc, ^bc} + min{ZAc, ^bc} > L, 
and we use the auxiliary events 

£' = {2Zac + 2bc > L}, 
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and 

to bound ¥[8] as follows 

P[^'] < P[^] < U £"] < 2P[^'], 

where we used that ¥[8'] = ¥[£"] (by symmetry) in a union bound, and the fact 
that, on 

2max{ZAc, 2bc} + minlZAc, 2bc} > 22ac + ^bc > L. 

Hence 

--lnP[£'] > --lnP[£] > --ln2P[£'] = --lnP[£'] - - In 2, 
L L L L L 

and, taking a limit as L — )■ +00, 

aR.(t)= lim -^lnP[^]= lim -^lnP[f'], 
provided the hmit exists. 



Moment-generating function In order to compute the limit above, we use the 
moment-generating function 

0(.)=E[exp(s[24'^ + 4')])], 



(which does not depend on i) as described in Section 2.3 Dividing up the expec- 
tation into the four possible cases, we have 

= (^p + ^{1 - p)^ + ^{1 - p){e' + e^') < +00, 

for all s e M. Letting 

W, = ^il-p), 



the derivative of (h(s) is 
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Decay rate By large-deviations theory, we are looking for a solution to 

Letting = e^, we get the quadratic equation 

(p + Wp) + Wp{uj + cu^) = Wp{uj + 2uj^), 

or, rearranging, 

{p + Wp)^Wpiu^ 

whose solution is 



Ip + Wp 

Wp ■ 



Then 



a^*{t) = lim -^lnP[£:'] 

Z/->+oo L 



lim lnP[2ZAc + ^bc > L] 
s* — ln0(s*). 



Noting that 



2{p + W,) + ^{p + Wr)W^, 



we get 

s*-ln0(s=^) = In 
Rearranging, we have finally 



^/p + Wp 



'^^/Wp{p + Wp) + w^p VF+w^ 7 ' 



-ln|2W^e-(l-^e-)+^e- 
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Asymptotics By a Taylor expansion, we get as i — > that 



an*{t)^^t' + 0{t'). 



On the other hand, as i — > +00, 



-In e-*/2 2 



1 f 1 _ He-.) + le-/^ 
3 V 3 y 3 



where 




A.3 STEAC/SC 



Definitions For a locus i, we let D)^^ be the time to the most recent common 
ancestor of A and B in Ge (in units of N generations). We let 



Similarly, we define D'll., D^'^, Vac and Pbc- Then STEAC/SC fails if 

Dab > min{PAc,I^Bc}, 

an event we denote by £. Once again, to deal with the term on the RHS, we 
re-write S as 

Dab - min{PAc, I^bc} > 0, 
and we use the auxiliary events 

S' = {Vab - Vac > 0}, 



L 



Vab = J2^- 



and 



S" = {Vab 



Vbc > 0}, 



to bound F[S] as follows 



F[S'] < F[S] < F[S' U S"] < 2F[S'], 
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where we used that P[£^'] = P[£^"] (by symmetry) in a union bound, and the fact 
that £' impUes S. Hence 

-\\n¥[£'] > -j\nF[S] > -jln2F[£'] = -ylnP[£'] - yln2, 

Ij Ij Ij Ij Ij 

and, taking a limit as L ^ +00, 

Q;sTEAc(i)= lim -ylnP[^]= lim -^lnP[£:'], 
provided the limit exists. 



Moment-generating function In order to compute the Umit above, we need the 
moment-generating function 

0(.)=E[exp(.[Dj) -Di^)])], 

(which does not depend on £). Dividing up the expectation into the four possible 
cases, we have 

+ ^(l-p)E[e-^^] 
+ ^(l-p)E[e^^^] 
4(1 -P) 



where we used: 



1. In the case SUCCESS^, D^^g — tab = Eq where Eq is an exponential mean 1 
conditioned to be less than t. Independently, using the memoryless property 
of the e 
Hence 



of the exponential, D\q — tabc = Eq where £'0 is an exponential mean 1. 



^AB - ^AC = -^AB + ^0 - TABC - E^ ^ -t + - Eq. 

2. In the case Fail^ and T[Gi] = AB|C, D^H^ - tabc = ^1 where Ei the 
minimum of (2) independent exponentials mean 1, that is, an exponential 
mean l/Q = 1/3. Moreover, ofl — Tabc = Ei + Ei where Ei is an 
exponential mean 1 independent of £^1. Hence 

^AB - ^AC = ^ABC + ^1 - TABC - E, - E, = -E,. 
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3. In the case Fail^ and r[G^] = AC|B, of^ — tabc — El where Ei is an 
exponential mean 

l/@ = 1/3. Moreover, of^ — tabc — Ei-\- E-i where 
El is an exponential mean 1 independent of E-i. 

DaI - = TABC + Ei+Ei- TABC - E^ = E^. 

4. In the case Fail^ and r[Gi] = BC|A, D^l = D^H. 
Note that 

1 5 

for all |s| < 1, and 

1 . l-e-(i-^)* 



E[e^-^°] = - /" e^^e-^dx = 
P Jo 



Hence 



(s) = pe 



P Jo P(l - s) 

1 _ e-(i--)* 1 



-St ■ 



p{l — s) 1 + S 



1 - S2 3 V 1 - S2 



3(1 - s2) ■ 

The derivative of (f){s) is 

_ [-3fr-'-2.sr-q[3(l-,s^)]-[3r-'-A-q[-6s] 



-t 



[18s - 9t{l - s^)]e-'^ - 6se 
[3(1 - s^)]^ 

Decay rate By large-deviations theory, we are looking for a solution to 

0'(s) _ [6s - 3t(l - s2)]e-"* - 2se-* 
(l-s2)(3e--*-s2e-*) ' 
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Note that the denominator on the RHS is positive on s G (0,1), and that 

m 

and 

lim = +00, 

s-^l- (p[S) 

SO that by [|Dur96ll there is a solution < s* < 1 to ([3]). The solution must 
satisfy 

[6s, - 3t(l - sl)]e~''' - 2s,e-* = 0. (4) 
which can be re- written as the fixed-point equation 

s, = ^[6s,-3til-sl)]e^^-''^' = Ft{s,), < < 1. (5) 

Note that Ft{0) = -3te* < and Ft{l) = 3 > 1. Moreover, 

Fl{s) = ^[6 + 6ts]e(i-^)*- ^[6s-3t(l-s')]e(^-^)* 
= ^e(^-^)*[6 + 3^2(1 - s^)] > 1, 

for < s < 1. Hence Ft{s) — s is strictly increasing and has a unique solution in 
(0, 1). Eq. ([5]) is easily solved numerically. 
Then 

asTEAc(i) = lim -^lnP[,f:'] 

= lim -ylnP[pAB-2^AC >0] 
= — ln0(s*). 

Asymptotics We consider asymptotic s when t — )■ 0. Define 
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Evaluating the LHS in Q (for e > small but fixed) as t — )■ gives 

[6s^ - 3t(l - s^)]e~"=* - 2s,e-^ 
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9 



t + Qe-'r -3t{ 1 r + 0(r) 

\ 16 



l-t + - + 0(f) 



so that, because 



-1 



< and As 



-1 



> 0, 



the solution of (|4]) satisfies s_e < < for < £ < | and t small enough. 
Then 



3e — s^e * 



3(1 - si) 

3 (i - + o{e 



9 

16' 



X 
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1 + r 



■Mi 



Since this holds for all e > small we get 



e + o{e)] [i-t + - + o{e) 



asTEAc(t) = -ln0(s*) = -e + 0{e). 



For the t — )► +oo asymptotics, let 

1 



■u = - and cr = (1 — 



23 



Substituting in (|5]), we get 



1-(TU = ^[6(1 - (Tu) - 3a(2 - au)]e'', 
which after rearranging becomes 

Se" -1- Sae" 

u 



Sae"" - a - ^a'^e'^ 



(6) 



We have J-'(O) = +oo. Moreover, letting cr* be the only positive solution to 

g{a,) = 3e"* - 1 - 3a,e"* = 0, (7) 

we have J'(cT,) = 0. Note that ^' (a) = -3ae'^ < 0,^(0) = 2 andlim^^+oo ^(a) = 
— oo, so that cr* is well-defined. Noticing that Q appears in the denominator of ([6]) 
as well we get that is strictly decreasing between a = and a = a*. Hence the 
limit t — 7- +00 is equivalent to the limit a ^ . Finally, in that limit, letting at 
be the a-value giving rise to the value u = 1/t 

asTEAc(i) = -ln0(s*) 

3g-..t _ ^2 -t 



-In 



t - In 



t - In 



3(1 - si) 
3e"* -il-^Y 



t 



3(l-(l-f)2) 



t-\nt-/3t, 



where 

lim = In = In — = cr* — In 2, 

t->+oo 6cr* 6cr* 

where we used (|7]). 
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